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Abstract. Point defects such as interstitials, vacancies, and impurities in 
otherwise perfect crystals induce complex displacement fields that are of long- 
range nature. In the present paper we study numerically the response of a two- 
dimensional colloidal crystal on a triangular lattice to the introduction of an 
interstitial particle. While far from the defect position the resulting displacement 
field is accurately described by linear elasticity theory, lattice effects dominate in 
the vicinity of the defect. In comparing the results of particle based simulations 
with continuum theory, it is crucial to employ corresponding boundary conditions 
in both cases. For the periodic boundary condition used here, the equations of 
elasticity theory can be solved in a consistent way with the technique of Ewald 
summation familiar from the electrostatics of periodically replicated systems of 
charges and dipoles. Very good agreement of the displacement fields calculated in 
this way with those determined in particle simulations is observed for distances 
of more than about 10 lattice constants. Closer to the interstitial, strongly 
anisotropic displacement fields with exponential behavior can occur for certain 
defect configurations. Here we rationalize this behavior with a simple bead-spring 
that relates the exponential decay constant to the clastic constants of the crystal. 
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1. Introduction 

The properties of crystalline substances often crucially depend on the structure and 
dynamics of imperfections of the crystal lattice. In particular, point defects such as 
interstitials and vacancies play a pivotal role in determining the stability, transport 
properties, growth characteristics, and mechanical behavior of materials. Recent 
impressive experimental advances, such as optical tweezers and confocal microscopy 
[H [2] , now permit to study the fundamental properties of point defects in condensed 
matter systems with "atomistic" space and time resolution. 

Recently, a number of experimental studies have focused on the structure and 
dynamics of point defects in two-dimensional assemblies of micrometer sized colloidal 
particles [3l|4l[5] and, in particular, on their effective interactions [6l[7l[8]. In studying 
such defect interactions the question arises to which degree they can be rationalized in 
terms of continuum elastic theory. As a first step towards answering this question, in 
this article we investigate numerically the disturbances caused by isolated interstitial 
particles and compare the results with the predictions of continuum theory. In 
carrying out such a comparison, it proves crucial that in solving the equations of 
elasticity theory boundary conditions are used that match those of the simulations. 
For the periodic boundary conditions usually applied in computer simulations, the 
displacement fields of single defects can be determined using the technique of Ewald 
summation familiar from electrostatics [HI [TO]- While elasticity theory properly 
describes the disturbances and interactions created by lattice imperfections on a larger 
scale, discrete lattice effects dominate on spatial scales of the order of few lattice 
constants. 

The remainder of this paper is organizes as follows. In Sec. [2] we define the 
model and describe the numerical methods. The treatment of point defects in a 
two-dimensional elastic continuum is discussed in Sec. [3] and comparison with the 
numerical results is discussed in Sec. ID For certain defect configurations one observes 
an exponential rather than algebraic decay of the displacement fields. This behavior 
can be understood in terms of a simple bead-spring model introduced in Sec. [5] with 
parameters related to the elastic constants of the material. Some concluding remarks 
are provided in Sec. [51 

2. Simulations 

In this paper we study a two-dimensional crystal of soft particle interacting via the 
Gaussian potential [TTl [TH [H] 

u(r) = eexp(-r^/cr^), (1) 

where r is the inter-particlc distance and e and a set the energy and length scales, 
respectively. In the following, energies are measured in units of e and distances in units 
of a. This so-called Gaussian core model, used here as a generic model for a system 
of soft spheres, is a realistic description for the short-ranged effective interactions 
between polymer coils in solution |14| . In three dimensions, the Gaussian core model 
can exist as a fluid, a bcc- and an fcc-solid depending on temperature and density 
[12j . In two dimensions, the perfect triangular lattice is the lowest energy structure of 
Gaussian core particles at all densities |15] . Computer simulations indicate that also in 
this system of purely repulsive particles point defects such as interstitials, vacancies or 
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impurity particles of different size display attractive (as well as repulsive) interactions 
both in two and three dimensions [16j . 

To study the displacement field of a single interstitial numerically, we prepare a 
configuration of particles arranged on the sites of a perfect lattice configuration and 
insert an extra particle of the same species. After insertion, the system is relaxed to 
a new minimum energy configuration by steepest descent minimization, i.e., we study 
the defect structure at T = 0. Typically, 70.000 steepest descent steps are carried out. 
The system we study here consists of = 199.680 Gaussian core particles (without the 
extra particle) at a number density of p = 0.6a~^ corresponding to a lattice constant 
a = (2/V3/9)^/^ ~ 1.3872(7. Periodic boundary conditions apply to the simulation box 
of length = 416a and height Ly = (^3/2)480 a = 415.692a. The aspect ratio of 
the almost square simulation box is Ly/Lx = 0.99926. 

We quantify the perturbation caused by the defect in terms of the displacement 
field [m 

u(rO EE r: - r,. (2) 

Here, and denote the position of particle i with and without the defect, 
respectively. As we will see in the following sections, simple point defects generate 
remarkably intricate displacement patterns that can be understood in terms of 
elasticity theory only on large length scales. 

At T = 0, the elastic constants describing the macroscopic response of the system 
to perturbations can be calculated as a function of density from simple lattice sums. 
For a density oi p = Q.Ga^, the Lame coefficients (see Sec. [3]) of the perfect triangular 
lattice have values A = 1.1487 ecr~^ and ^ = 0.06018 ecr"^. At this density, the pressure 
is p = 0.5442 ea and the energy density is e = 0.2691 ecr~^ corresponding to an energy 
per particle of E/N — 0.4485e. The bulk modulus, which in two dimensions is related 
to the the Lame coefficients by if = A + p, has a value of = 1.2089ecr~^. 



3. Elasticity Theory 

While close to a point defect the displacement field is highly anisotropic and strongly 
dependent on the atomistic details of the interactions, for large distances elasticity 
theory is expected to be valid. The differential equations describing the equilibrium 
of an elastic continuum are usually expressed in terms of the strain tensor [17j 

1 / dui du-j \ 

^-(^)-2(^^^)' 
where Ui denotes the i-component of the displacement u and the i-th component of 
the position r. For a given external volume force f(r) with components ft acting on 
an isotropic system such as a crystal on a triangular lattice. Hook's law leads to the 
equilibrium condition for the strain: 

Here, A and /i are the so-called Lame coefficients and summation over repeated indices 
is implied. Solving this equation for a singular force yields the Green's function from 
which the response of the elastic continuum to an arbitrary force can be obtained by 
integration. 

To model the displacement field caused by the introduction of point defects using 
linear continuum elasticity theory, we determine the displacement field caused by two 
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pairs of opposing forces, one pair acting along the x-axis and the other one along the 
y-axis [TFl [IH HO] ■ This idealized model of a defect is equivalent to inserting a small 
circular inclusion into a hole of different size |19| . Each force of this pair is of equal 
magnitude F but with opposite sign acting on two points separated by a small distance 
h. Such a force pair exerts a zero net force on the material. In the limit /i — > where 
the force _F — > oo in a way such that Fh remains constant, the equilibrium condition 
for the displacement can be written as 

{X + 2^,)^ = Fh6ir). (5) 

Assuming that the displacement can be written as the derivative of a potential. 



one obtains 



A(-^0)=-^(r). (7) 



This equation is the Poisson equation of electrostatics with a singular disturbance. 
Since, as noted above, K{i) = — ln(r)/27r is a solution of AK = — <5(r) (see, for 
instance, Rcf. [21]), we obtain the Green's function 

from which the displacement field u(r) follows by differentiation according to Equ. 

m, 

" 27r(A + 2^) r2 ' 

In comparing the results of particle simulations with those of elasticity theory it is 
important to realize that the displacement fields predicted by continuum theory are of 
a long-range nature. Therefore, it is crucial that corresponding boundary conditions 
are used in both cases. All simulations discussed in this paper are done with periodic 
boundary conditions in order to minimize surface effects and preserve the translational 
invariance of the perfect lattice. Hence, also the continuum calculations need to be 
carried out with periodic boundary conditions. 

Since the defect fields for the infinite material are long-ranged, the displacement 
field in the periodic system cannot be obtained by simply summing up the 
contributions of the periodic images. In fact, such a naive summation of the 
contribution of all image defects diverges. A more appropriate treatment that avoids 
this problem consists in determining the Green's function of the Poisson equation 
([7]) for periodic boundary conditions. In this case, the solution of this equation in 
two dimensions, known from electrostatics [9l I10|. can be written as Ewald sum of a 
logarithmic potential embedded in a neutralizing background. 



27r(A + 2/i) 



2^^e-^'/4''' ^, ^ TT I 
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Here, Ei{x) = Jl^{e^/t) dt is the exponential integral. The first sum is over all lattice 
vectors 1 in real space and the second sum is over all reciprocal vectors k in Fourier 
space. The adjustable parameter 77, set to a value of ry = Q/L^ here, determines the 
rate of convergence of the two sums and A is the area of the rectangular simulation 
cell. The Fourier space sum can be evaluated accurately using about 2,500 reciprocal 
space vectors. From Equ. (fTO|) for the scalar function (j){r) the displacement field of a 
point defect in a system with periodic boundary conditions is found by differentiation. 



For the systems considered in this paper, the real space sum may be truncated after 
the first term. Since the value of F/i/27r(A + 2/i) is undetermined, the parameter 
7 = Fh/2iT{X + 2fi) is treated as a fit parameter in the following. The Ewald sums 
of the above equations describe the effects of "image defects" at the center of the 
periodically replicated domains. 

4. Results 

First, we study the displacement field of a single interstitial. To generate such a 
defect, we insert an extra particle of the same species into a perfect 2d-crystal on a 
triangular lattice. After insertion, the system is relaxed to a new minimum energy 
configuration by steepest descent minimization, i.e., we study the defect structure at 
T = 0. Typically, 70.000 steepest descent steps are carried out. In each step each 
particle is moved in the direction of the force acting on the particle where the absolute 
value of the displacement in chosen to be small enough to ensure that the energy of 
the system decreases in each step. 

The extra particle can deform the crystal in different ways [H and produces 
displacement fields of different symmetries (see Fig. [l}. In one configuration, called 
I2 interstitial or crowdion and shown in Fig. [T^, the additional particle pushes one 
particular particle of the crystal out of its equilibrium position. Both the original 
particle and the additional particle arrange themselves at equal distance around the 
lattice position of the original particle. The displacement pattern arising for this type 
of interstitial has two-fold symmetry and, of course, occurs in all three low-index lattice 
directions with equal probability. One may suspect that this defect configuration, 
with a symmetry that differs from the symmetry of the underlying triangular lattice, 
is caused by the rectangular periodic boundary conditions that are applied to the 
system. To rule out this possibility, we have repeated the calculation with hexagonal 
periodic boundary conditions obtaining the same result. 

Another low-energy defect configuration is the interstitial with three-fold 
symmetry (sec Fig. [TId). In this case the interstitial particle is located at the center of 
a basic lattice triangle and pushes its neighbors outward from their original positions. 
A third important interstitial configuration is the Id interstitial or dumbbell interstitial 
shown in Fig. [TJ:. In the 2d Gaussian-core model under the conditions studied here 
the I2 pattern has a slightly lower energy than the I3 interstitial and the Id interstitial. 
The energy difference between a I2 and a I3 interstitial is 0.000674e and the difference 
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Figure 1. Displacement fields (top) and local defect configurations (bottom) for 
the I2 (a), the /a (b) and the l^i (c) configurations. The length of the arrows 
representing the displacements of the particles from their position in the perfect 
lattice are exaggerated for better visibility. In the figures at the bottom the blue 
spheres represent the particles and the small gray spheres indicate the position of 
the lattice sites of the perfect crystal. 



between I2 and Id is 0.000665e. All three displacement patterns are important for 
the diffusion of interstitials. An I1 interstitial is very mobile in the direction of its 
main axis. The /a and forms are visited as intermediate configurations when the I1 
interstitial changes the orientation of its main axis and hence its direction of motion 
[22]. 

Next, we compare the displacement fields determined numerically with the 
predictions of elasticity theory. In particular, we verify to which extent the 1/?'- 
behavior modulated by the periodic boundary conditions and embodied in Equ. (|lip 
is realized in the particle system. The complex displacement patterns of the various 
interstitial configurations shown in Fig. [1] obviously differ from this expectation, at 
least near the defect, and indicate that continuum theory is not applicable in this 
region. Far away from the defect, however, the perturbation caused by the defect 
is small and the response of the material should be described accurately by linear 
elasticity theory. 

The magnitude |u(r)| of the displacement vector u(r) is shown as a function of the 
distance from the interstitial in Fig. [2] for the li defect configuration. Each point in 
the figure corresponds to one individual particle. For short distances, the displacement 
magnitude is not a unique function of the distance r reflecting the anisotropic nature 
of the defect. For larger distances, however, the displacement magnitude is mostly 
determined by the distance r. Eventually, however, the periodic boundary conditions 
lead to a spread of the displacement magnitude for even larger distances and a 
splitting into two branches corresponding to the x- and y-directions and the directions 
along the diagonals, respectively. In the regime where u(r) behaves isotropically, the 
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Figure 2. Displacement magnitude |u(r)| as a function of distance from the 
defect r for the I2 interstitial. Each red dot corresponds to one particle. The 
solid line represents the 7/r behavior. Here, 7 = 0.229lcr^ was used as this value 
yields the best fit of the results obtained vie Ewald summation to the results of 
the particle simulations in the far field. Inset: angle 9 between the displacement 
vector u the position vectors r as a function of the distance r from the defect site. 

displacement follows the approximately 1/r-form predicted by elasticity theory for a 
point defect in an infinitely extended medium. The orientation of the displacement 
vector u(r), depicted in the inset of Fig. [21 behaves in an analogous way. The angle 
between u(r) and the position vector r, shown as a function of the distance r from 
the defect, is not a unique function of r near the defect. For larger r, 9 vanishes 
indicating that in this distance regime the displacement vector points straight away 
from the defect. At even larger distances, the periodic boundary conditions imposed 
on the system eventually cause the angle 6 to spread again. 

The displacement fields calculated according to Equ. pT|) and numerically for 
an interstitial in the I2 configuration are compared in FiglS) In this figure, the 
displacement components Ux and Uy are depicted as a function of the distance from the 
defect along the x-axis and y-axis, respectively. The prediction of continuum theory, 
calculated using the Ewald summation of Equ. PT|) . agrees well with the displacement 
field of the particle system for distances larger than about 10 lattice constants. 

5. Harmonic Model 

Near the defect, the predictions of continuum theory differ from the simulation results. 
The deviation is particularly pronounced in the direction of the main axis of distortion 
of the I2 defect, in which the displacement appears to decay exponentially up to a 
distance of about « 10 a. This unexpected exponential behavior can be understood 
in terms of a simple model with harmonic interactions. This model consists of a one- 
dimensional chain of particles in which each particle is connected to its two neighbors 
with springs of force constant fci (except the first and last particle, which are coupled 
only to their neighbors on the right and left, respectively). In addition, each particle 
is attached to a fixed lattice position with another spring of force constant k2. The 
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Figure 3. Displacement components Ux and Uy of the I2 interstitial as s function 
the of the distance along the x-axis and y-axis, respectively (solid lines). Here, 
the direction of largest displacement of the I2 defect is oriented in x-direction. 
Also plotted is the displacement computed from continuum theory according to 
Equ. mil l (dashed line), simple 1/r-behavior (dotted line) and the displacement 
obtained for the simple mechanical model described in the main text (dash-dotted 
line). The inset shows the region close to the defect location. As in Fig. [2]a defect 
strength of 7 = 0.2291tT^ was used for the evaluation of the displacement from 
elasticity theory. 



Hamiltonian of this system is 

7^ = I - - 5)^ + f E - ^^■)'' (12) 

where + 2 is the number of particles, Xj is the position of particle j and h is 
the equilibrium distance between two neighboring particles. In the minimum energy 
configuration of this chain, the particles arc arranged such that Xj = jb. We now 
imagine that particle is pushed to the right by a distance of uq while particle + 1 
is kept fixed at xn+i = {N + 1)6. If the system is then relaxed to a new energy 
minimum, all other particles will be displaced from their original positions too. For 
this simple model, the response of the system to the displacement of the first particle 
can be calculated analytically by direct matrix inversion (see [Appendix A[ ) . In the 
large N limit, one finds that the displacement of the particles from their original 
position decays exponentially with their position, 

Uj = Moexp(-Q;j), (13) 

where Uj is the displacement of particle j due to the forced displacement uq of the 
first particle. The decay constant a is related to the force constants of the model by 

a = cosh^i (^1 + . (14) 

To compare the prediction of this simple model with the simulation results we 
have to determine the force constants ki and ^2 felt by the particles in the main 
axis of the defect. While the force constant ki arises from interactions within this 
main axis, the force constant k2 is related to interactions of the particles in the main 
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axis with those from adjacent rows. Accordingly, we determine ki by calculating 
numerically the energy change caused by slightly displacing one single particle in a 
one-dimensional row of otherwise fixed Gaussian core particles without the presence 
of the neighboring rows. The distance of the particles in the row is chosen to be equal 
to the lattice constant at the density p = 0.6a~^ considered throughout the paper. 
From the energy as a function of the displacement one obtains a force constant of 
ki = O.OlSe/fj^. To determine the force constant k2 we calculate the energy change 
caused by translating a whole row of particles in the perfect crystal. The particles in 
the row are fixed with respect to each other and the remaining particles are kept at 
their lattice positions. From the energy change per moved particle a force constant of 
k2 = 0.0013e/(j^ follows. The decay constant of a = 0.29 calculated according to Equ. 
(|14p with these force constants is in perfect agreement with the computer simulation 
results shown in Fig. [3l 

For a system in which only nearest neighbor interactions are important, the force 
constants ki and ^2 can be simply related to the bulk modulus K and the shear 
modulus /i. Then, the force constant ki is given by 

fci = 2v"{a), (15) 

where v(a) is the pair potential at distance a. Since in this case the elastic moduli are 
given by 

(16) 

and 



2 ' ' a 



2 / S/i , 



one obtains 

To the extent that the response of the system to shear is determined by the interaction 
of neighboring parallel rows of particles, the energy density caused by shifting a whole 
row of atoms between two fixed ones is the same as that of a shear of appropriate 
magnitude. Accordingly, the force constant k2 is related to the shear modulus by 

k2 = ^/i. (19) 

This expression remains also valid if interactions beyond nearest neighbors are included 
between adjacent rows of particles. In terms of the elastic constants, the constant a 
describing the exponential decay of the displacement field along the principal axis can 
be expressed as 

or, in terms of the Poisson ratio 

a = cosh-(l^). (21) 

For a density of p = 0.6(t~^, inserting the Poisson ratio of v = 0.905151 determined 
from a simple lattice sum yields a « 0.25, only slightly different from the correct 
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value a w 0.29. This deviation occurs, because in the Gaussian core model at the 
density p = Q.Q(j~^ interactions between non-nearest neighbor particles are important 
in determining the elastic constants (in fact, considering only nearest neighbors would 
produce a negative shearing modulus /i in this case). For systems, in which only 
nearest neighbor interactions are relevant, the above expression is expected to hold 
accurately. 

6. Conclusion 

Point defects in two-dimensional crystals, such as interstitials and vacancies, can 
assume configurations with symmetries that vary from the symmetry of the underlying 
triangular lattice. While close to the defect the displacement field is highly anisotropic 
and strongly dependent on the atomistic details of the interactions, for large distances 
elasticity theory, which predicts isotropic behavior, is valid. For the particular 
I2 interstitial configuration, the displacement decreases exponentially with distance 
along the main defect axis. The decay constant is simply related to the material 
properties via the Poisson ratio, which measures the ratio between transversal and 
axial strain upon stretching. In comparing the displacement fields computed from 
particles simulations with those obtained with continuum elasticity theory it is crucial 
to use equivalent boundary conditions in both cases. Since particle simulations are 
usually carried out with periodic boundary conditions, also the differential equations 
of elasticity theory need to be solved for a periodic system. We have shown here that 
Ewald summation, a technique routinely used in computer simulations to determine 
the electrostatic interactions of charges and dipoles, can be used for this purpose. In 
this method the sum over all interactions with periodic image defects is split into 
two sums in real space and reciprocal space, respectively. This particular treatment 
of the long-ranged nature of displacement fields effectively introduces a neutralizing 
background that leads to convergent sums. Note that exactly the same expression 
apply also to a system that is enclosed in a rigid container. Outside a core region 
near the defect, displacement patterns determined using such Ewald summation agree 
perfectly with those calculated in particle simulations. 
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Appendix A. 



nix) 



In this appendix wc calculate the response at T of the one-dimensional bead- 
spring model of Sec. [5] to a forced displacement of the first particle in the chain. The 
potential energy of the TV -I- 2 particles, located at positions Xj, is given by 

N , N+1 

i=o j=o 
where b is the equilibrium distance and fci and k2 are force constants. The vector 
X = xq, xi, ■ ■ ■ , xn+1 includes the positions of all particles. Minimizing the potential 
energy with respect to the particle positions xj by requiring that 

dH{x) 



fci 
2 



dx-i 



= 



(A.2) 



for all j, one finds that at the potential energy minimum the particle positions are 
Xj = bj. We now displace particle by an amount uq from its original position afo = 
and keep particle -I- 1 fixed at position (N + 1)6. If we hold particle at this new 
position while minimizing the potential energy, all particles from 1 to A will move to 
new equilibrium positions. Thus, the minimum energy configuration of the system is a 
function of the displacement ito of particle 0, which may be viewed as a parameter that 
is controlled externally and perturbs the system. To make this distinction between 
the displacement of particle and that of all other particles more explicit, we denote 
Uq with an extra symbol, ^ = uq. The displacement Uj of the particles j ~ 1, • • • , A 
is then a function of ^, 



u,iO^Xj{0-x,iO), 



(A.3) 



where xj (^) 



and Xj{0) denote the particle position in the minimum energy 



configuration with and without perturbation, respectively. In the following, we will 
calculate the displacements Uj(^) as a function of the perturbation strength ^. 

Since condition (jA.2p defines the position of the energy minimum as a function 
of the perturbation strength ^, its derivative with respect to ^ must vanish. 



0. 



d_ ( mm 

Application of the chain rule then leads to 

d^n \ dxm 



E 



dxj dxi 



This condition must hold for all i. Defining 



dxjdxi 



x=x{^) 



and 



d^dxi 



(A.4) 
(A.5) 

(A.6) 

(A.7) 
(A.8) 

(A.9) 



x=x(£) 



Displacement fields of point defects in two-dimensional colloidal crystals 



13 



we can rewrite Equ. (jA.Sp as 

y, = Y,n,,z,. (A. 10) 

j 

Inversion of the matrix Titj then yields the vector z, 

j 

Once Zj = cfxj{£,)/d£, is known, Xj{£,) can be obtained by integration. 

For the bead-spring model considered here, the first and second derivatives of the 
potential energy with respect to the particle coordinates are given by 

{2ki + k2)xi- kiXi+i - kiXi^i - k2bi, (A. 12) 

if j = J, 

if j = i + 1 or j = ?: - 1, (A. 13) 




^^i^^J I else 



and 



(A.14) 



du _ r -fci if i = 1, 

dxidi ^ \ if i > 0. 

To solve Equ. (jA.lip we have to invert the symmetric tridiagonal matrix Hij. For 
this particular matrix, the inverse matrix is known analytically |23j . 

1 fcosh[(A^ + l-b--i|)«] 



^ ~ 2ki\ sinh(a) sinh[(iV + l)aj 
cosh[(iV +1-1- j)a] 
sinh(a) sinh[(A^ + l)a] 



(A.15) 



where 



a = cosh-i (l + ;^ ) (A.16) 



(A.18) 



2fci 

Since for our model y — {— fci, 0, 0, • • • , 0}, we obtain 

^-T.HIi'yj=-H^,'k,. (A.17) 
j 

and hence 

dxi _ cosh[(7V + 2 - i)a] - cosh[(iV - i)a] 
~ 2sinh(a) sinh[(iV + l)a] ' 

For large N , this equation simplifies to 

^=exp(-*a). (A.19) 

Integration with respect to ^ then yields 

5,(0 =eexp(-ia) + a, (A.20) 

where the integration constant is given by Ci — Xi{Q). Thus, the displacement of 
particle j is proportional to the displacement of particle and decays exponentially 
with the distance from the origin, 

Ui — uoexp(— za), (A. 21) 

with a decay constant a that depends on the force constants ki and fc2 only. 



